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Abstract. We incorporate a time- independent gravita- 
tional field into the BGK scheme for numerical hy- 
drodynamics. In the BGK scheme the gas evolves via 
an approximation to the collisional Boltzmann equa- 
tion, namely the Bhatnagar-Gross-Krook (BGK) equa- 
tion. Time-dependent hydrodynamical fluxes are com- 
puted from local solutions of the BGK equation. By ac- 
counting for particle collisions, the fundamental mecha- 
nism for generating dissipation in gas flow, a scheme based 
on the BGK equation gives solutions to the Navier-Stokes 
equations: the fluxes carry both advective and dissipa- 
tive terms. We perform numerical experiments in both ID 
Cartesian geometries and axisymmetric cylindrical coor- 
dinates. 

Key words: hydrodynamics — methods:numerical — 
shock waves — gravitation 



1. Introduction 

The BGK code distinguishes itself from other hydrocodes 
in that it has recourse to the physics which generates dis- 
sipation, namely the physics of particle collisions. Devel- 
oped by Prendergast & Xu (1993), the hydrodynamical 
scheme invokes a microscopic description of gas flow and 
it is therefore based on considerations of gas kinetic theory. 
Recall that kinetic-based hydrocodes rely on the fact that 
the state of a gas can be described by giving the distribu- 
tion function of particle velocities /(x, u, t) at a point in 
the phase space of a single particle. The codes take advan- 
tage of the fact that the quantities of hydrodynamical in- 
terest (namely the mass, momentum and energy densities 
in the gas) are low-order velocity moments of /. / obeys 

D f t) f co ^ 

the Boltzmann-equation, which we write as -gt = -gr , 
where is a time rate of change along the trajectory 
of a single particle moving freely in phase space (under 
the action of smoothly varying forces, if these are present) 

and is the rate of change of / due to collisions. 

The classical two-body collision integral is non-linear, and 
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non-local in velocities (Cercignani 1988). BGK (from the 
names of Bhatnagar-Gross-Krook) replaces this integral 
with the term ^ g ~^ , where g is the Maxwell-Boltzmann 
distribution function having the same mass, momentum 
and energy densities as /, and r is a relaxation time, which 
can (ideally) be as small as the mean time between colli- 
sions of a particle in the gas (Bhatnagar, Gross & Krook 
1954). It is important to notice that / and g have the 
same mass, momentum and energy densities, but they do 
not give rise to the same fluxes of these quantities. Fluxes 
are determined by different moments of the distribution 
functions. BGK differs from all other kinetic-based codes 
in several respects: the collisions are active throughout the 
duration of a time step, and are not imposed as a distinct 
"equilibration" process at the end of a timestep. Also, we 
do not guess the form of the distribution function but 
solve the BGK equation to find it. The BGK equation is 
linear in /, which makes it easy to solve if g is known. 
However g is not known; the "compatability" conditions 
that / and g have the same mass, momentum and energy 
densities would determine g if / were known. Therefore 
the BGK formulation is also nonlinear and non-local in 
velocity space, as is the Boltzmann equation. It will be 
shown that this situation leads to a set of non-linear inte- 
gral equations for the parameters of g. It might seem that 
we have made no progress by replacing the Boltzmann 
equation by the BGK plus compatability conditions; but 
this is not true, because it will be shown that we need 
solve for the parameters of g only in the neighborhood of 
a boundary between computational cells, and for a short 
time (given by the usual CFL (Courant-Friedrichs-Lewy) 
condition). Knowledge of the parameters of g is equivalent 
to knowledge of the mass, energy and momentum densi- 
ties. 

The connection between the BGK equation's micro- 
scopic description of gas flow and a macroscopic descrip- 
tion has been shown by Cercignani (1988) for the case of 
a perfect monatomic gas and by Xu (1993) for polyatomic 
gases. Velocity moments of the BGK equation give the 
Euler equations for negligible particle collision time, r, 
the Navier-Stokes equation for small yet non-zero r and 
a description of rarefied gas dynamics for large r. In the 
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Navier-Stokes regime, these derivations also furnish ex- 
pressions for the shear and bulk viscosity coefficients and 
the heat conductivity coefficient. 

In principle, if the local value of the collision time can 
be measured as a function of time for a gas, then the BGK 
scheme may be used to evolve the gas with true physical 
dissipation parameters. In practice, the resolution of the 
BGK scheme is limited by the grid resolution. If the grid 
is not fine enough to resolve a discontinuity, then artificial 
dissipation must be added to broaden the discontinuity 
so that it is at least one grid cell thick. Because viscosity 
and heat conductivity are proportional to r, the BGK 
scheme broadens shocks by enlarging r at the location of 
the discontinuities. The expression for the collision time in 
the BGK scheme therefore contains two terms. One term 
is the real physical mean collision time and it is chosen 
according to the desired Reynolds number of the problem. 
The second term is chosen in such a way that shocks in 
the flow span at least one grid cell. The latter term tunes 
the amount of artificial dissipation in the scheme. The 
notable difference between how the BGK scheme inputs 
artificial dissipation and how other schemes input artificial 
dissipation is that the BGK scheme puts it in exactly as if 
it were real dissipation corresponding to the numerically 
necessary value for r. 

We emphasize that modelling the real dissipation in 
an astrophysical object such as that resulting from "tur- 
bulent" viscosity, is still outside the reach of any existing 
hydrocode. For one thing, values for viscosities in astro- 
physical objects are highly speculative. Secondly, to solve 
Navier-Stokes problems with a particular value of the vis- 
cosity for an astrophysical object requires grid resolutions 
which are beyond what is currently achievable. For lab- 
oratory scale phenomena and for small enough Reynolds 
number, the BGK code is successful at modelling the dis- 
sipation produced by particle collisions in a fluid. Tests of 
the Kolmogorov and the laminar boundary layer problem 
show real dissipative effects (Xu & Prendergast 1994). 

Thus far the BGK code has been extensively tested 
without the inclusion of gravity. Several papers document 
the tests which verify its accuracy and robustness as both 
an Euler and Navier-Stokes solver in ID and 2D Cartesian 
geometries and for two-dimensional adaptive unstructured 
grids. The list of one-dimensional Euler tests which have 
been performed with the BGK scheme includes: the Roe, 
Sod, Lax-Harten, Woodward-Colella, and Sjogreen tests 
for subsonic and supersonic expansion (Prendergast & Xu 
1993; Xu 1993; Xu, Martinclli & Jameson 1995). The re- 
sults from these one-dimensional test cases are that the 
BGK scheme produces shock fronts which are typically 
one to two cells wide, and contact discontinuities which are 
slightly broader. This is competitive with high resolution 
codes which do not employ regridding in the neighborhood 
of a shock front. The BGK scheme also exhibits negligi- 
ble under and over shooting even at strong shock fronts. 
The other notable features in the tests with rarefaction 



waves are the sharp corners at the junctions between the 
rarefaction waves and the undisturbed, uniform regions. 

Many codes have trouble treating a low density re- 
gion with a flow gradient. The BGK scheme successfully 
handles low density regions because it satisfies both an en- 
tropy condition and a positivity condition (Xu et al. 1996). 
In the BGK scheme collisions relax the gas toward local 
thermodynamic equilibrium states and this relaxation pro- 
cess is accompanied by an increase in entropy. Godunov- 
type schemes on the other hand typically demand an en- 
tropy fix when they encounter strong rarefaction waves (cf. 
the Sjogreen test) otherwise they produce unphysical rar- 
efaction shocks. Usually the fix is an addition of artificial 
viscosity. Because the BGK scheme naturally satisfies the 
entropy condition, it simply cannot generate these unphys- 
ical phenomena. In satisfying the positivity condition, the 
BGK scheme avoids producing states with negative den- 
sity or internal energy. The Roe test is conducive to the 
creation of these non-physical states but the BGK scheme 
does not encounter them (Prendergast & Xu 1993). The 
same cannot be said of the performance of conservative 
finite difference schemes and codes employing Riemann 
solvers (e.g. Roe's approximate Riemann solver). 

The list of two-dimensional Euler test cases which have 
been performed with the BGK scheme includes: (a) uni- 
form Mach 3 flow in a tunnel with a forward facing step 
(the Emery test). The BGK code is not modified in any 
way near the step to treat the flow past it and its corner, 
which is a singular point. With the BGK scheme expan- 
sion shocks never emerge from the corner (Prendergast & 
Xu f 993; Xu 1993: Xu et al. 1995), (b) double Mach reflec- 
tion in supersonic flow over a wedge, (c) the diffraction of 
a strong shock (Mach = 5.09) around a corner. Test cases 
(b) and (c) are further examples of the BGK scheme suc- 
ceeding at simulating flow without having to summon any 
detection algorithms or entropy fixes. According to Xu et 
al. 1995, the original Godunov scheme, the Roe scheme 
without the entropy fix, and the Osher scheme could pro- 
duce a rarefaction shock at the corner in test case (c) ; (d) 
flow around an impulsively started cylinder. When applied 
to this problem, many schemes either fail or have severe 
problems in maintaining positive pressure and density in 
the near vacuum low pressure and low density region cre- 
ated behind the cylinder. The BGK scheme seems to be 
able to preserve positivity without ad hoc fixes, and to 
reach a steady state solution for the problem (Xu et al. 
1996). 

To test its performance as a Navier-Stokes solver and 
to see if it has real viscosity effects, the BGK scheme has 
been applied to the laminar boundary layer problem and 
the Kolmogorov problem (Xu 1993). The laminar bound- 
ary layer problem models the flow of gas above a flat plate. 
Even on coarse grids (e.g. (32 X 16), (16 X 8)), the BGK 
scheme impressively recovers the Blasius profile. In the 
Kolmogorov problem, a one-dimensional sinusoidal veloc- 
ity field is imposed in a uniform density and isothermal 



A. Slyz & K.H. Prendergast: BGK with Gravity 



3 



fluid. The BGK code fulfills the Navier-Stokes prediction 
which is that the shape of the fluid's velocity profile re- 
mains unchanged while the amplitude of the velocity de- 
creases in such a way that the fluid's kinetic energy decays 
exponentially. The agreement between the theoretical vis- 
cosity coefficient and the numerical viscosity coefficient 
(deduced from the measured decay rate) is excellent. 

To improve the resolution of physical discontinuities 
occurring in complicated flows, a version of the BGK 
scheme on a two-dimensional adaptive unstructured grid 
has now been developed (Kim & Jameson 1998). 

Most astronomical applications of a hydrocode require 
a consideration of gravity. To this end we have been de- 
veloping the BGK scheme. Recently the BGK scheme has 
been used for a cosmological simulation (Xu 1997) but this 
is prior to results showing the long-term stability of the 
BGK scheme with gravity and the BGK scheme's conver- 
gence to the equilibrium state with gravity. In this paper 
we give such results. In addition to the incorporation of 
gravity, we have modified the geometry of the Eulerian 
grid to axisymmetric cylindrical coordinates. Changing 
the geometry of the grid has effects similar to the effects of 
adding gravity to the BGK scheme. It gives rise to source 
terms in the hydrodynamic equations. For the purpose of 
simplicity, in this paper we describe modifications to the 
BGK method when a time-independent gravitational po- 
tential is incorporated into a BGK scheme in Cartesian 
coordinates (section ||). We present results however for 
both gas flow on a one-dimensional Cartesian grid and on 
an axisymmetric cylindrical grid (section ||). 

2. A BGK Flow Solver with Gravity 

2.1. Hydrodynamic Equations from the BGK Equation 
with Gravity 

In Cartesian coordinates, we specify the position and ve- 
locity of a particle with coordinates (x, u). To be clear, x 
here is a 3- vector (x, y, z) giving the position of a particle 
in configuration space and u is a 3-vector (u,v, w) giving 
the velocity of a particle in this space. In this coordinate 
system and in the presence of a smooth gravitational field 
with gravitational potential, $(x), the BGK equation is: 



df 8f 8$ df _g-f 

dt <9x 9x du t 



(1) 



(For reference, this is equivalent to equation 2.15 in Shu 



(1992) where % is the BGK collision term 2=£.) Re- 
call that in the BGK collision term g is the equilibrium 
distribution towards which the true velocity distribution 
function / relaxes on a collision timescale, r. We define 
the units of / and g to be the mass in a volume element 
of phase space. The particle collision time r depends on 
macroscopic quantities such as temperature and density, 
and it is therefore a function of position x and time t. 
Both / and g are functions of space x, time t, and par- 



ticle velocities u and £ where £ = £j> • • • £k) is a K- 
dimcnsional vector of velocities associated with the inter- 
nal degrees of freedom of a particle. With this parame- 
terization of the internal velocities, the energy associated 
with internal degrees of freedom is written as \m^ 2 where 
£ 2 = £i 2 + £ 2 2 + • • • + S.K 2 - Included in these internal de- 
grees of freedom is the energy not explicitly accounted for 
in versions of the code with physical spatial dimension, 
D, less than three. For example a particle in a diatomic 
gas (7 = 7/5) moving in 3 spatial dimensions has 5 de- 
grees of freedom. (The connection between the number of 
degrees of freedom, n, and the ratio of specific heats, 7, 
is 7 = (n + 2)/n.) Because a one dimensional hydrocode 
(D = 1) treats only one of these degrees of freedom ex- 
plicitly, the behavior of the diatomic gas is reproduced by 
treating the four remaining degrees of freedom as internal 
energies. We represent those energies in the K-dimensional 
vector £ where K = n — D. In terms of 7, K = ^zrq — D. 

Equations governing the time evolution of the mass 
density, momentum density and kinetic plus internal en- 
ergy density are found by multiplying (|l|) in turn by 1, u, 
and 5 (u 2 + £ 2 ) and integrating over a volume in velocity 
space, gE = du^ K ~^ d£. The volume element in veloc- 
ity space has the above form in the £ variable because 
we need only the total internal energy of a particle and 
not the energy associated with any particular internal de- 
gree of freedom. Hence an integration has been performed 
over angles in £. An essential point about the derivation 
of the hydrodynamic equations from velocity moments of 
the BGK equation is that moments of the BGK collision 

term, i2—Jl ) vanish. This is enforced by the fundamental 
hypothesis of the BGK model: in the relaxation process 
collisions tend to reshuffle the particles defining the true 
distribution function, /(x, u, t) into a local equilibrium 
distribution, g(x, u, t) having the same mass, momentum 
and energy densities as /(x, u, t). We call this requirement 
for the equivalence of / and g's moments, the "compata- 
bility condition" . 

For detailed results of taking velocity moments of (|l]), 
we refer once again to Shu 1992 (p. 20-23). We give the 
moment equations here in compact form: 



d_ 

dt 



(0 



d 

dxk 



dx k duk' 







(2) 



where we have defined the following bracket notation: 
(...)= J(...)fd3 

and C = 1 , , (u 2 + £ 2 ). In eq. |^ we see that mass, 
momentum and energy densities are updated by terms 
which are the divergence of a quantity (i.e. this quantity 
is commonly referred to as a flux) and when gravity (or 
a curved coordinate system) is involved, also by terms 
which cannot generally be written as the divergence of a 
quantity and hence cannot update the contents of a cell 
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by a flux through the cell's surface. The latter terms are 
called source terms. Without gravity and in Cartesian co- 
ordinates, a BGK scheme for hydrodynamics evolves all 
the hydrodynamic quantities exclusively by flux terms. We 
emphasize that there is no dissipation source term in equa- 
tion H Dissipative terms are carried by the fluxes. The 
BGK scheme is unique in this respect. In Cartesian coor- 
dinates it allows a fluid to evolve concurrently through ad- 
vective and dissipative processes without decoupling them 
into two separate operations. 

With gravity, source term computation is unavoidable. 
The gravitational source term in the momentum equation 
cannot be manipulated into the form of the divergence 
of some quantity, so it persists as a source term. The 
gravitational source term in the energy equation on the 
other hand may be reformulated to give an energy equa- 
tion without a source term if the gravitational field is in- 
dependent of time. In this paper we consider this case. 
Without a gravitational source term, the energy equation 
for the case of a time-independent gravitational field is a 
conservative equation for the total kinetic, internal and 
gravitational energies. In one dimension it is: 

Q Q 

(Skin + Stat + p$) + {F Skb 



* — &v- ^W)=0. 

We define F £ki " +£int here as the flux of the kinetic plus 
internal energy, and F p as the mass flux. In F £ki » +£int 
and F p we are only considering the x-component of the 
flux, i.e. the flux through a wall perpendicular to the x 
direction. 

In the end, the BGK scheme follows the time evolution 
of the integrated values of mass, momentum and energy 
densities within cells. For example, for one dimensional 
flow in Cartesian coordinates the BGK scheme tracks the 
contents of a volume element with boundaries x\ and X2- 
If the values of mass, momentum and energy densities (de- 
noted below by q) are given within this cell at time t\ then 
at time ti the new values of the mass, momentum and en- 



ergies are updated by fluxes {F) (section 2.2) through the 
cell's bo und aries and in some cases by source terms (S) 
(section 2.7). 



X2 pX 2 

q(x,t2)dx= / q(x,ti)dx + 

X\ J X\ 

*2 /"*2 

T(x\ , t) ■ x dt — / J-(x2,t)-xdt + 



S(x, t) dx dt. 



2.2. Hydrodynamical Fluxes with Gravity 

Fluxes in the BGK scheme arise from velocity moments 
of a particle distribution function / which is a solution to 
the ordinary differential equation + - = - with initial 
conditions / — f at t — t = 0. This equation holds along 
each trajectory. The solution is composed of two terms: an 



integral over the past history of g and a term representing 
relaxation from an initial state, /q. 



/(x/,u,t) = - f 3 (x,u',i')e-( t - i ')/^<' + 
T Jo 

e~ t/T / ( X/ - uf,u,i ) 



(4) 



Here x and u in the arguments of g are solutions of a 
gas particle's equations of motion in Cartesian coordinates 
with gravity: 



dx 
~dY 



du d$ 

dt' dx 



with the final conditions x = Xf and u = u at t = t. To 
first order in time the solutions of these equations are: 



u(t — t ); u 



t). 



(5) 



These equations define a single trajectory for the parti- 
cles which arrive at Xf at time t with velocity u. The 
trajectory would generally be curved if we carried terms 
of higher order than those shown. Even when we are for- 
mulating the BGK scheme for a I-dimcnsional Cartesian 
grid there is no escape from considering u a vector in 
three dimensions. Both / and g always depend on all of 
the "molecular" velocity components because individual 
particles can have velocities of any speed and in any di- 
rection. A consequence of this is that when fluxes are even- 
tually computed at xj, /(x/,u,i) is integrated over the 
full and continuous range of velocities, u, from — oo to 
+oo. Therefore fluxes computed from the true distribu- 
tion function /(x/,u, t) arise by a weighted integration of 
the equilibrium distribution function g(x , u , t ) not just 
over one trajectory but over all possible trajectories (in 
6-dimensional pha se s pace) which arrive at (xf,t). (We 
leave it to section 2.6 for a discussion of the implications 
of this.) For a I-dimensional problem in Cartesian coor- 
dinates all we require is that the trajectory wind up on 
the wall located at (xf,yf,Zf) where yf and Zf are fixed 
for all x because the only non-trivial dimension is the x- 
dimcnsion. 

Now as already stated in the introduction, g is not 
known a priori along all the trajectories passing through 
x/ at time t. Physically g must describe an equilibrium 
(3) state so it is sensible to assume a Maxwellian form for g: 



g(x',u.',t')=p 



, x {K+3)/2 

M -A((u'-U) 2 +C 2 ) 



(6) 



Here U = (U, V, W) are the mean macroscopic veloci- 
ties in the x, y and z directions respectively. To explain 
why U has this physical interpretation we go to the one- 
dimensional problem. In one-dimension, V and W are zero 
so g is effectively left with 3 parameters: p, U and A. For 
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a perfect gas with K = 2/(7 — 1) — D the "compatability" 
condition gives an explicit relationship in closed form be- 
tween the mass, momentum and energy densities and g's 
parameters: 




9 

pU 
(C/ 2 + ^ 



-3) 



(7) 



2 A 



(Note that when g is Maxwellian its velocity moments fol- 
low simple recurrence relations (see appendix in Xu et al. 
(1996)).) Eq. (j] assigns physical meaning to g's parameters: 
p is the mass density, U is the mean macroscopic velocity 
in the x direction, and A is related to the mass density, p, 
and the internal energy density, £ m t, hi the following way: 



A = 



(K + 3) p 



where K was defined earlier as the dimension of the vector 
£. A is therefore inversely proportional to the temperature 
of the gas. Of course mass densities, mean macroscopic 
velocities and internal energy densities change with space 
and time in the fluid and their evolution is assumed to be 
governed by the evolving distribution function /. Hence 
there is a mutual connection between / and g. By substi- 
tuting the formal solution for / (eq. [|) which itself depends 
on g into the "compatability condition" on the moments 
of g, one arrives at a set of coupled non-linear integral 
equations for the space and time dependence of g's pa- 
rameters: 



^ ae -X^ a (r,t) dE = 



10 



T f ipafo{x-f - ut,t )dE. 



(8) 



Here we have adopted the following abbreviated notation 
for g and its parameters. Let 

g = e-x^° 

where a = 1 , . . . , 5 and 
a = (l,u',i(u' 2 + e 2 )). 

We define the transpose of ijj a to be -0 Q T and \a to be a 
vector of g's parameters: 

Xa = Hn[p(X/TT) (K+3 ^%-2XU,-2\V,-2XW,2X). 

Notice that ip a has no space dependence whereas Xa has 
no velocity dependence. In equation (||) we assume that 
the collision time, r, which is a function of macroscopic 
densities is locally constant and hence can be removed out- 
side the integral which is performed over the velocit y mo- 



ments. The computation of r is discussed in section |2.3.2 . 
Because fo in the formal solution for / (eq.^) may be 



straighforwardly constructed from the initial conditions, 
and because the integral over g(x', u',t') in the formal so- 
lution for / is readily evaluated once g(x', u', t') is known, 
the bulk of the numerical work in the BGK scheme lies 
in solving the above coupled non-linear integral equations 
for 5 (x',u',t'). 

2.3. Numerical Computation of f 

Following the sketch of what is required for a computation 
of /, we outline some aspects of the numerical procedure, 
particularly those which are relevant to adding gravity to 
the scheme. The details of the flux computation for the 
2D Cartesian case without gravity are presented in Xu 
& Prendergast (1994), Xu et al. (1996), Kim & Jameson 
(1998). In Cartesian coordinates we divide the computa- 
tional domain into cells. We then suppose that we are 
given the values of mass, momentum and energy densities 
within each cell at the beginning of a timestep, to. We 
describe the procedure for computing / at a point on 
a cell boundary. 

2.3.1. Numerical Computation of g 

We begin the numerical computation of the true distri- 
bution function /(x^,u, t) at x/ by the construction of 
<?(x , u , t ) which appears in the integrand of (Q). Because 
we only need to know / at position x/ located on a cell 
interface and for the duration of a CFL timestep we do 
not compute g(x , u , t ) precisely at each point along all 
the trajectories crossing xy at t. Instead we compute g 
at the point x/ on the cell interface at time to (where to 
corresponds to the beginning of a time step) and then ap- 
proximate it in the spatial and temporal vicinity of (xj, to) 
through a Taylor expansion. 



g(x',u',t') = 5(x/,u,t )(l 



dhig 



(x'-x/) 



dhig 
Ijt 7 ' 



(f - to) + . . .) 



(9) 



To show the velocity dependence in the spatial and time 
derivatives of g we rewrite the Taylor series expansion as: 



ff(x ,u ,t ) =g(x / ,u,i )(l + a- (x -x f ) 

+A.(t'-t )). 



(10) 



Here the vectors a and A depend on the derivatives of 
the parameters of g with respect to space and time and 
on u and £ 2 as well. The spatial and time derivatives of 
g are assumed to be locally constant and we group them 
as coefficients of particle velocities, u , £ in the following 
way: 



/ / / / ,2 

a = ai + a-iu + CL3V + 041V + as I u + £ 

A = Ax + A 2 u + A 3 v + A 4 w' + A 5 (u' 2 + £ 



(11) 
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The subscripts 1, 2, 3, 4, and 5 refer to the five parameters 
in g (cf. eq. ||). Recall that x and u lie on curved trajecto- 
ries (eq. ||) if the trajectory is sufficiently prolonged. How- 
ever if we substitute the expressions for x , u given in (||) 
into the expansion for g(x , u , t ) (eq.nG) (with the expres- 



sions for a and A given in (11)) and if we retain only terms 
linear in time, the curvature terms in the particle trajec- 
tories do not contribute to the expansion. Hence to first 
order in time, the solution of the BGK equation for the 
true distribution function /(x/,u,£) in one-dimensional 
Cartesian coordinates with gravity is equivalent to the so- 
lution for the one dimensional BGK equation in Cartesian 
coordinates and without gravity. This result generalizes to 
the case of non-Cartesian geometries. It is the first impor- 
tant result concerning the incorporation of gravity into 
the BGK scheme. 

Our strategy for solving for the terms in the Taylor 
expansion of g is: 

(a) solve for g(x.f, u, to) and the first spatial derivatives of 
g from initial conditions. This is possible because as we 
already showed in section 2.2 there is a direct connection 
between g's parameters and the mass, momentum and en- 
ergy densities which are specified within each cell at the 
beginning of a timestep t . We can easily obtain the values 
of the macroscopic densities and their first spatial deriva- 
tives with respect to x' at (x/,io) through interpolation 
and thereby solve for g and g's first spatial derivatives 
at the cell boundary, x/. We refer to Prendergast & Xu 
(1993) for the details of the numerical interpolation in- 
cluding a discussion of interpolation switches. The latter 
are devices to prevent spurious maxima and other arte- 
facts, and are commonly used in the reconstruction phase 
in many schemes. 

(b) before solving for the time derivatives of g, we con- 
struct the initial distribution function, fo- Because the 
distribution function / should reflect a possible non- 
equilibrium initial state across the cell interface located 
at position x/, we assume fo to be composed of two half 
Maxwellians. It is easiest to see what we mean by this by 
again considering a one-dimensional Cartesian problem. If 
we divide the computational domain into intervals whose 
boundaries are specified by xi,X2, . . ., x n and we look 
at flux transfers in the x-direction through one of these 
boundaries located at Xf then we may construct fo at Xf 
by constructing one half-Maxwellian to the left of xj and 
another half-Maxwellian to the right of Xf. The param- 
eters of the two half Maxwellians as well as their slopes 
are obtained from the initial conditions, just as g's pa- 
rameters were obtained but the left (right) Maxwellian is 
computed from the mass, momentum and energy densities 
interpolated from the left (right) side of the cell interface. 
This choice for the form of fo is not unique. The details of 
fo should not matter, since they are rapidly damped for 
small r. 

(c) Finally we solve for the time derivatives of g by in- 
sisting that the "compatability condition" is satisfied on 



average over the CFL time interval. Due to our approx- 
imate computation of g and hence of /, g and / cannot 
have exactly the same moments for the duration of a CFL 
time step over all the trajectories arriving at x/. However 
we insist that they have the same moments at x^ and on 
average over the time interval < t < t for which the 
true distribution function is computed. 



J J ^a T (/(x/,u,f') -g(x f ,u,t'f) d^dt 







(12) 



This condition gives an expression for A (see Xu & Pren- 
dergast (1994)). It is worth remarking that the manner 
in which the time derivatives of g are computed is in 



some sense implicit. Because eq. 12 is applied over the 
entire time interval, expressions for A arise from infor- 
mation throughout the entire time interval as opposed 
to just information at the beginning of a time step. We 
also point out that by finding the time-dependence of / 
through the formal solution for / plus the "compatabil- 
ity condition" applied on average over an updating time 
step, Prendergast and Xu's (1993) BGK scheme bypasses 
the stiffness in the BGK equation which would force the 
updating timestep to be the particle collision timescale 
r. Because the updating time step is determined from the 
CLF condition hydrodynamics may be performed with the 
BGK scheme. Earlier attempts at implementing the BGK 
equation (e.g. Chu 1965) were constrained by the stiffness 
of the BGK equation and their application was therefore 
limited to rarefied gas dynamics. 

2.3.2. The Collision Time, r 

Aside from the Courant factor in the computation of the 
CFL time step (section |2.5| ) (and certain constants which 
depend on the adopted interpolation rules), the expression 
for the collision time, r, contains the only 2 parameters 
in the code. According to gas kinetic theory, the mean 
time between collisions in a gas with density p and tem- 
perature T is proportional to l/(pT 1 / 2 ). Since A w l/T, 
t = C\\f\j p where C\ is a proportionality constant and 
the first of the two parameters in r. We add a second term 
to the expression for the collision time. In one-dimensional 
Cartesian coordinates this second term is: 



Pi - \TKI Pr\ \Pl~Pr\ 



M/ pi + VK/p r {Pi +Pr) 



(13) 



The subscript I (r) denotes quantities interpolated from 
the left (right) side of the cell interface at which we com- 
pute the collision time r and p is the gas pressure. We 
choose the above form for the second term in the expres- 
sion for the collision time r because it helps with shocks 
across which there are gradients in p. As stated in the in- 
troduction (section [l]), we motivate this second term in the 
collision time as follows: given that the coefficients of heat 
conduction and viscosity are proportional to the collision 
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time t, the second term in the expression for r acts to in- 
crease both heat conduction and viscosity at shock fronts 
but not at contact discontinuities across which pressure 
is continuous. When the second term in the collision time 
dominates, the order of the scheme is reduced and the true 
distribution function / is determined more by the initial 
distribution function / than by the integral over g. 

2. 4- Computing the Total Fluxes 

Once f(x.f,t) is fully known as a function of the parame- 
ters of g, fo, the collision time r, and t, the time-dependent 
fluxes at Xf can be computed. We illustrate flux computa- 
tion for a one-dimensional Cartesian problem. In that case 
we are interested only in the x-component of the total flux, 
i.e. the component of the flux which is perpendicular to a 
cell wall at Xf. The total amount of mass, x-momentum 
and energy densities transferred through the wall perpen- 
dicular to the x-direction in a CFL time step < t < T is 
found from: 

TP 

^r£kin + £int 



/ f I 1 \ 

u x f(x f ,u,t)dEdt (14) 

Note that in the end, the gas kinetic origin of the expres- 
sions for the fluxes disappears and they are written as 
highly non-linear functions of the mass, momenta and en- 
ergy densities which are presumed known at the beginning 
of the time step. 

2.5. Courant Time Step 

The time step T used to evolve the dynamical equations 
is chosen so that it satisfies the CFL condition. Physi- 
cally the CFL condition limits the distance that informa- 
tion can travel in one time step to one cell. For the one- 
dimensional Cartesian case we compute T for each cell on 
the grid from 

Oi+l/2 T 2 + 
(\Ui+l/2\ + Cj+1/2) T + £ i+ i/2 = 

where a is the local gravitational acceleration due to grav- 
ity, i.e. — ^jr, U is the local macroscopic velocity, c is the 
local adiabatic sound speed and 

A+1/2 = min((x l+ i - Xi), (x t - X;_i), (xj_i - x 4 _ 2 )). 

This choice for Ci+1/2 insures that flux transport between 
adjacent cells is carried out for a length of time short 
enough so as not to empty them. 

Finally the time step T is chosen from 

T = (1 - e)min(T) (15) 

where (1 — e) is a safety factor; it is called the CFL factor 
(typically ess.6) . 



2.6. Comments about Flux Computation with the BGK 
method 

Unlike Ricmann solvers which propagate information 
along the characteristics of the Euler equations, the BGK 
scheme propagates information along the characteristics 
of the Boltzmann equation. These are particle trajecto- 
ries in phase space. In the BGK scheme fluid properties 
along the full continuum of trajectories passing through 
(x, t) contribute to the instantaneous fluxes at (x,t). The 
implications of this are many. Firstly it bypasses the limi- 
tation of lattice Boltzmann codes which discretize velocity 
space thereby restricting the number of directions for the 
propagation of information and thereby also limiting the 
magnitudes of the gas velocities . 

Secondly, using all the trajectories with appropriate 
weighting endows the BGK scheme with intrinsic upwind- 
edness. An upwind scheme computes fluxes using infor- 
mation coming from the same direction as the flow. When 
computing fluxes, the BGK scheme does not select trajec- 
tories corresponding to the direction of flow and discount 
the rest; instead it uses information carried by all trajec- 
tories arriving at the place where we construct /, but it 
weights the information along each trajectory according to 
the number of particles assigned to the trajectory by the 
BGK solution. This is what we mean by an intrinsically 
upwind scheme. 

Thirdly, the consideration of the full continuum of tra- 
jectories for the flux computation creates the potential for 
a truly multidimensional code. Because the velocity in- 
tegration includes all propagation directions and because 
the trajectories can be followed multidimcnsionally, the 
flux in one direction can potentially be computed from in- 
formation coming with all speeds and all directions. For 
a genuinely multidimensional scheme, the interpolation of 
the initial conditions for the reconstruction stage must 
also be multidimensional. For example, for a computation 
on a two dimensional Cartesian grid, the reconstruction 
of the initial conditions at a point (x, y) should include 
information from cells surrounding the point in both the 
x and the y directions. 

At the present time, the multidimensional BGK 
schemes which have been implemented do not take ad- 
vantage of the truly multidimensional possibility that the 
BGK scheme offers. They still simplify flux updates by 
separating them into individual operations in independent 
directions. This is an example of "operator-splitting" . On 
the optimistic side, a truly multidimensional BGK scheme 
is conceivable. The same cannot be said for schemes based 
on Ricmann solvers. Even in two dimensions it is not feasi- 
ble to construct a truly multidimensional Riemann solver. 
Upon representing the initial conditions as piecewise con- 
stant in two-dimensional cells, the Riemann solution is 
straightforward to apply at each of the interfaces of a two- 
dimensional cell but problematic at cell corners. Further- 
more, by propagating information only along the direc- 
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tions normal to cell interfaces, a two-dimensional Riemann 
solver has the same limitation as an operator-split method: 
it does not accommodate the possibility that waves in 
two-dimensional flows may propagate in infinitely many 
directions (Roe f986). 

2. 7. The Source Terms 

Because to first order in time the computation of the true 
distribution function, f(x,u,t), and hence the computa- 
tion of the hydrodynamical fluxes is unaffected by gravity, 
gravity's influence on the gas flow is relegated to source 
terms. 

2.7.1 . Gravitational Momentum Source Terms 

In one-dimensional Cartesian coordinates gravitational 
source terms contribute to the momentum change in a 
cell over a CFL time step t = to t = T in the following 
form: 



S Px (x i ^ 1 /2,T) = - 



p(x',t')^-(x',t')dx'dt'. (16) 
ox 



This integral may be computed in a number of ways. One 
possible approach is to do a linear interpolation in x' using 
the values of the mass densities and forces at the cell wells 
surrounding the cell for which we wish to compute the 
source term. 

The time integration required by ( |l6| ) takes advantage 
of the fact that we know the new mass densities at the 
end of an updating time step, i.e. at t = T, before we've 
completed the updates of the other fluid quantities, i.e. P Xl 
£ kin + £'mt + £grav ■ This is a consequence of there being no 
source terms in the mass continuity equation. From the 
mass densities at times t = and t = T we may perform a 
bilinear interpolation in x and t for the integral in equation 
M). Specifically 



Sp x { x %-x/2,T) = 

T r Xi 9$ 

■- J {p{x i - 1 ,0) — (x i - 1 ,0)(l-^)) + 

p(x i ,0) — (x i ,0)(Kx')) + 

pix^T) — (x^ u T)(l-l3(x')) + 
d<S> 



p(xi,T) — (xi,T)(p(x'))}dx' 



(17) 



where 



P(x') = 



(x' - Xj) 

(x i+1 -Xi)' 



2.7.2. The Energy Equation 

As for gravity's contribution to the energy equation, we 
explore two versions of it. One form of the energy equation: 



~bt 



(£kin + £int) 



d_ 

dx 



P: 



9$ 

dx 







-P x §f . The gravi- 



has a gravitational energy source term, —P x ^. 
tational energy source term may be computed analogously 
to the way in which the gravitational momentum source 
terms are computed (cf. eq. 17). 

As mentioned in section 2.1, another form of the en- 



ergy equation incorporates the gravitational energy source 
term into the energy fluxes and hence is in conservative 
form: 



dt (£kin 



d 

■ £int + £grav) + (F k '" + 



<PF P ) = 0. 



The conservative form of the energy equation requires a 
good computation of the gravitational energy because now 
it, in addition to the kinetic energy, £kin, must be sub- 
tracted from the total energy for the purpose of obtaining 
the internal energy, £; nt . These subtractions are poten- 
tially extremely delicate. They are important because the 
internal energy is used in computing A (which is inversely 
proportional to the temperature), a crucial parameter for 
the hydrodynamic fluxes. 

Because the form of the gravitational energy is similar 
to the form of the gravitational source terms, 



p(x )&(x )dx 



(18) 



it may be numerically computed in the same manner, i.e. 
by a linear interpolation to the values of the densities and 
the potentials at the walls bounding a cell. 

For the gravitational flux term, ^(<&i ?p ), in the con- 
servative energy equation we use the value of the grav- 
itational potential $ at the same point Xf at which we 
construct the true distribution function, /, for the other 
fluxes. 



2.8. Summary of the BGK Method with Gravity 

We have given a description of the procedure for com- 
puting hydrodynamical fluxes and source terms in a BGK 
scheme designed for hydrodynamics in the presence of a 
time-independent gravitational field. 
1.) While incorporating gravity into the BGK scheme we 
have found that to first order in time gravity does not alter 
the computation of the distribution function from which 
hydrodynamical fluxes are computed. It enters into the 
hydrodynamical computation through gravitational mo- 
mentum source terms as well as through either a gravita- 
tional energy source term or through gravitational energy 
fluxes. The latter depends on the chosen form of the en- 
ergy equation. 
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2. ) The energy source term formulation by definition does 
not guarantee conservation of total (kinetic + internal + 
gravitational) energy. An inaccurate computation of the 
energy source term introduces numerical heating and/or 
cooling. Because the energy equation can be manipulated 
into a conservative form for the total kinetic + internal 
+ gravitational energy, computation of an energy source 
term is avoidable. In lieu of an energy source term, the 
gravitational energy enters into the hydrodynamical com- 
putation through flux terms. These flux terms do not re- 
quire any flux computation in addition to the computa- 
tion which is already required in a hydrodynamics scheme 
without gravity. This is because they arc conveniently the 
product of the gravitational potential and the mass flux 
terms. 

3. )Even though an implementation of the conservative 
form of the energy equation by definition assures to- 
tal energy conservation, it introduces another challenge 
brought on by the inclusion of gravity into a hydrocode. 
For the calculation of hydrodynamical fluxes in the BGK 
scheme it is essential to accurately compute the tempera- 
ture and therefore the internal energy. With gravity and 
with an implementation of the conservative form of the en- 
ergy equation it therefore becomes necessary to accurately 
compute the gravitational energy so that it, in addition 
to the kinetic energy, may be subtracted from the total 
energy (kinetic + internal + gravitational) to give the in- 
ternal energy without introducing numerical heating (or 
cooling) into the scheme. 

3. Results 

We present results from two tests of the BGK scheme 
in the presence of a time-independent gravitational po- 
tential. The first test reveals the effect on a simulation 
of using different forms of the energy equation. The sec- 
ond test is performed on an axisymmetric cylindrical grid. 
Both tests reveal the capability of the BGK scheme with 
gravity to reach an equilibrium state and to maintain a 
gas configuration in hydrostatic equilibrium. 

3.1. On a ID Cartesian Grid: Gas Falling into a Fixed 
External Potential 

We perform this test case to compare the results from 
two versions of the BGK scheme with gravity: one using 
the conservative form of the energy equation (the ECS 
scheme) and one using the energy source term formula- 
tion of the energy equation (the EST scheme). The initial 
conditions for the two simulations are identical. Each sim- 
ulation is performed on a ID Cartesian grid with 68 evenly 
spaced cells. The gas is initially stationary [P x = 0, where 
P x is the cc-momentum) and homogeneous with mass den- 
sity p = 1 and internal energy density £i n t = 1 in each 
cell. The external gravitational field is constant in time 
and its potential, <f> has the form of a sine wave, i.e. 



$ = -$o(£/(27r))sin(27ra/£) with $ = 0.02 and C = 64. 
Because we apply periodic boundary conditions there are 
two ghost cells at each end of the grid. We take e in the 
Courant condition to be .4, 7 = | and C\— .01 and £2= 1 
in the expression for the collision time r. We run both 
versions of the code for 500, 000 iterations. The total time 
at the end of the run with the ECS scheme is 285739.54 in 
machine units (4706 in units of the initial sound crossing 
time) . Note that we define the total time to be the sum of 
the lengths of each of the individual CFL timesteps used 
for the updates. For the run with the EST scheme the to- 
tal time at the end of the run is 261888.6 in machine units 
(4313 in units of the initial sound crossing time). 

Figs, [l] through ^ show results from the version of 
the BGK scheme which implements the energy conserv- 
ing form of the energy equation and which therefore keeps 
track of the total kinetic -I- internal + gravitational en- 
ergy density within each cell. Figs. ^ through [)] on the 
other hand show results from an implementation of the en- 
ergy equation in a non-conservative form, i.e. one in which 
gravity is incorporated into the energy equation through a 
gravitational source term. The plots show the time evolu- 
tion of various quantities (i.e. mass density, x-momentum 
and lambda) as functions of position on the grid. Because 
the cells are evenly spaced with a spacing of Ax — 1, we 
label cell positions by an index X which runs from 1 to 
68 in the plots. The quantity on the vertical axis in each 
of the plots is also plotted at approximately equal time 
intervals. 

For the simulations using the energy conserving 
scheme we present not only plots showing the long term 
evolution of the gas, but also some (Figs. |l|, |[ ||) which 
show the evolution of the gas early on in the simulation. 
In these figures, particularly in the x-momentum plot (i.e. 
Fig. ||), we see the steepening of the sound waves into 
shocks. With the values for the internal energy density 
and the mass density given in the initial conditions, the 
initial sound crossing time for the gas is 60.72 (in machine 
units). For the data presented in figures |l], || and ^ the total 
time elapsed is 201.19 (in machine units). This is equiv- 
alent to 3.3 sound crossing times and indeed we see that 
the gas has undergone a bit more than three oscillations 
in these figures. 

In the plots showing the long term evolution of the 
gas (Figs. H to ||) , at the time of the earliest curve (time 
w 2816) the gas has experienced about 46 sound cross- 
ing times (we compute this using the value for the initial 
sound crossing time). The mass density in both the ECS 
scheme and the EST scheme at t = 2816 is high at the 
location of the gravitational potential minumum and low 
at the location of the potential maximum indicating that 
the gas has fallen into the potential well. Because the gas 
should convert some of its gravitational energy into in- 
ternal energy as it falls into the potential well, we look 
to see if the gas at the density maximum is hotter. As 
expected, the temperature at t = 2816 is highest (hence 
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lambda is lowest) at the location of maximum density and 
the temperature is lowest (hence lambda is highest) at the 
location of minimum density. For reference, lambda at the 
beginning of the simulations is 0.75. 

In the results from the scheme using the non- 
conservative form of the energy equation, there is a mono- 
tonic decline in lambda (Fig. ^|) at increasing times. Since 
lambda is inversely proportional to the temperature, this 
decline corresponds to a temperature rise which we believe 
is due to "numerical heating" . This heating persists with 
time and therefore deters the gas from reaching thermo- 
dynamic equilibrium. The heating is reflected in the value 
of the total energy (kinetic + internal + gravitational) 
summed over the entire grid at the end of the run: it is 
~ 10% greater than the initial total energy. 

In contrast, the results from the run with the en- 
ergy conserving implementation of the energy equation 
show a gas which reaches thermodynamic equilibrium: the 
lambda plot (Fig. ||) is level at late times with a mean 
value of .77 (deviation ss 1.52D — 04) and remains at the 
same value throughout the entire grid for as long as we 
run the simulation. The difference between the initial and 
final total energies on the grid is 4.3451D — 13. In addition 
to thermodynamic equilibrium, the gas reaches mechan- 
ical equilibrium as verified by the ^-momentum plot: at 
the final time plotted P x is —1.161?— 20 in the mean with 
a deviation of 2.04Z? — 04. We believe that the noisiness in 
the final P x state is a consequence of the run being per- 
formed without the use of interpolation switches for the 
flux computation. The analytic solution for the equilib- 
rium mass density profile in a fixed gravitational potential 
is p{x) = p(0)e~( kr )* where p is the atomic weight, mu 
is the mass of a hydrogen atom, k is Boltzmann's constant, 
T is the temperature, $ is the gravitational potential. Be- 
cause the temperature should be constant throughout the 
grid at equilibrium, a plot of the numerical values of ln(p) 
vs. $ displays how closely the gas has settled to equilib- 
rium. We plot these quantities for the results from the ECS 
scheme in Fig. |l^. When we fit a line ln(p)^ = a+b<&i to the 
numerical results then a = -2.46484^-02, b = -15.7485 
and the variance of a is 2.47615-E — 05 while the variance 
of b is 1.7509OE - 03. 

3.2. On an Axisymmetric Grid: Gas Falling into a Fixed 
External Potential 

We simulate the shock heating of gas falling into a fixed, 
external gravitational potential well in 2D. The hydro- 
dynamic computation is performed on an axisymmetric 
cylindrical grid which has 50 logarithmically spaced cells 
in radius, w, and 100 evenly spaced cells in the z direction. 
The gas falls into a gravitational potential well (Fig. |l2|) 
which is derived from a spherical density distribution with 
the following profile: p ~ (1 + r a )~ 2,5 - The potential is 
centered on the origin of the axisymmetric cylindrical grid, 
namely at w — and at the midplane in z. The bound- 



ary conditions for the simulation are reflecting. The gas 
is initially stationary and uniform in density (Fig. |Tl| ). As 
seen in Fig. |ll|, it is initially higher in temperature at the 
location of the potential well than in other parts of the 
grid. In machine units, G = 7, p = 10, the radial mo- 
mentum and vertical momenta densities in each cell are 
0, i.e. P ro = and P z — 0. The internal energy density in 
each cell is initially set to £i„t = 2.5X10~ 14 — £ gra v, where 
£ g rav is the gravitational energy density in each cell. The 
axisymmetric cylindrical grid extends from —1.1 to +1.1 
in the vertical direction and from to r max = +1.1 in the 
w direction. We also take e in the Courant condition to 
be .8, 7 = §j and C\= .001 and €2= 1 in the expression 
for the collision time t. 

Figs. [TJI and |TJ show values of the density, the inverse 
of the temperature (lambda), and both radial and verti- 
cal velocity along two rays through the cylindrical grid. 
One ray goes through a point (vo = 0, <j) = 0, z = 0) and 
extends in the w direction. We plot various quantities as 
a function of w along this ray at time t n . The second ray 
goes through a point {m = zni, <j) = 0, z = —1.1), extends 
in the z direction and we plot various quantities as a func- 
tion of z along this ray at time t n . Thus in figs. [D| and [TJ 
each curve represents results at one instant (specified by 
t n ) in the course of a run. The curves are approximately 
evenly spaced in machine time units. The key for the lines 
styles can be found in the right hand corner of the density 
plot, (Fig. |l3[ ). In the following discussion we will refer pri- 
marily to these two figures ( |lj and |l^) and to the times 
for which information is presented in these figures. For a 
view of the gas evolution in the entire w — z plane, as 
opposed to merely along these two rays, we will also refer 
to a time sequence of density and lambda surface plots, as 
well as to plots of velocity vectors (Figs. E^, |l6, 17). In the 
surface plots, values of density and lambda are given on 
the vertical axis as a function of w and z in the horizontal 
plane. In the velocity plots, each arrow represents the vec- 
tor sum of the radial and vertical velocities at a point on 
the w — z plane. The length of the arrow is proportional 
to the strength of the velocity field. 

From ii through t-j, the gas density (Fig. |l^) rises at 
the centre of the grid, although the grid centre is not where 
it reaches its highest values initially (see also Fig. |l5|). This 
is because the strong gravitational forces (Fig. |l2|) keep the 
gas from falling uniformly into the gravitational potential 
well. Instead the gas density piles highest along a ring 
surrounding the grid's centre. The peak of the ring moves 
inward from ti through £7, and the density jump between 
the tip of the ring and the region outside of it is steepest 
and reaches a maximum at tg in these figures (see first 
density plot in Fig. |l^). Time tg marks the onset of the 
outgoing shock into the cold, low-density region outside 
the ring. 

The lambda plots (Fig. |TJ) show that between t\ 
and £7, the region outside the ring becomes progressively 
colder. This is a consequence of the regions nearer to the 
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ring falling inward more quickly then the regions further 
from the ring (Fig. |l4]) , leading to an expansion of the gas 
and a concurrent cooling. By i 7 the leading edge of the in- 
wardly spreading cold region has steepened into a shock. 
From ig through £15 the heated centre expands sending 
an outgoing shock which heats the gas as it propagates 
across the grid. We vividly see both the formation of the 
shock and its outward movement in the lambda surface 
plots (Figs. [l5|, [l6| and 17). In these surface plots, the 



cooling of the regions further from the ring is also visible. 

The formation of the shock is also clearly seen in the 
velocity profiles (Fig. [l4]) from t\ through t-j. During this 
time interval, the point of maximum velocity moves in- 
ward until at tg the ring material starts moving outward 
and encounters the material still infalling from the outer 
part of the grid at a shock front. From fg to tis we see 
the shock propagating outward. The outwardly moving 
shock front is delineated clearly in the velocity vector plots 
(Figs. and 0). 

After 100,000 iterations (Fig. |l|), the gas is almost 
in equilibrium: it has a smooth density profile and it is 
nearly isothermal and stationary. Computing the value of 

the free-fall time (tfj — yj 3 %Q p ) from the initial condi- 
tions, tff » .06. Therefore at 100,000 iterations (time 
= 36.3758) the time is 560 in units of the initial free-fall 
time. 

The most significant thing to notice in the results from 
this computation are the sharp shock fronts. Even though 
the computation was performed on a relatively coarse grid 
(50 cells in w and 100 cells in z), only about two cells are 
needed to resolve a shock and there is no evidence of under 
and overshooting. This is shown clearly in Fig. |l8|, where 
the inverse of the temperature, A, for tg — > ti2, is again 
plotted for the ray which cuts the cylindrical grid along 
the first cylindrical radius, w\ for fixed (j> = 0. However 
in Fig. we show the data for each of the four time 
instances in a separate plot, so as to better assess the 
shock resolution. The circles are the values of A in each 
cell and the thin lines are curves through these values. It is 
easy to see that only about two cells are needed to resolve 
a shock. 



4. Discussion 

The most notable results of the preceeding test runs may 
be summarized as follows: 

1. )For a fixed, external gravitational field both in Carte- 
sian geometry and in axisymmetric cylindrical geometry 
with a conservative form for the energy equation, the BGK 
scheme is able to keep a configuration in hydrostatic equi- 
librium for many sound crossing times. Furthermore with 
the energy equation in conservative form the BGK scheme 
settles to the correct physical solution. 

2. )Thc BGK scheme is able to achieve high resolution with 
rather coarse grids. For example the simulation of gas 



falling into a fixed external spherically symmetric grav- 
itational potential uses 50 cells in 137 and 100 cells in z 
and achieves shocks whose fronts span only 1 or 2 cells. 
We never resort to regridding in the neighborhood of a 
discontinuity to achieve high resolution. Conventional hy- 
drocodes frequently strive for high resolution through grid 
refinement techniques. This approach is both expensive 
and ineffective beyond a certain point. Uniform grid re- 
finement shortens time steps: in 3 dimensions, a decrease 
by a factor of 2 in the grid spacing of each of the 3 di- 
mensions involves a factor of 16 increase in the number 
of cycles that the code has to execute to reach the same 
total time. Furthermore even if one is willing to accept the 
expense of short timcstcps, no matter how fine the grid, 
a diffusive scheme never achieves the accuracy of a high 
resolution scheme run on the same sized grid. With that 
said, local adaptive grid refinement for the BGK scheme 
(without gravity) has been implemented by Kim & Jame- 
son (1998) and it gives excellent results for a number of 
tests involving unsteady supersonic flows. 

5. Conclusion 

We incorporated gravity into the BGK scheme for hy- 
drodynamics and presented results for the case of time- 
independent gravitational potentials both on a one- 
dimensional Cartesian grid and on an axisymmetric cylin- 
drical grid. Our results show that the BGK scheme is gen- 
eralizable to the case of a gas moving in the presence of an 
external, time-independent gravitational field. When we 
derive an approximate local solution to the BGK equa- 
tion, linear in space and time, gravity's curvature of par- 
ticle trajectories is unimportant. Gravity's effect on the 
flow enters into the computational method only through 
gravitational source terms. 

We conclude by stating that we believe that no other 
hydrodynamical scheme which is used in astrophysical 
computational fluid dynamics invokes as few "fixes" and 
operates with as much generality as the BGK scheme. 
There are no discontinuity detection algorithms, with sub- 
sequent special treatment of special regions. An entropy 
fix never has to be invoked when following the evolution 
of rarefaction waves because the BGK scheme satisfies 
the "entropy condition". Truly multi-dimensional BGK 
schemes (i.e. non-operator split) are realizable. 

We attribute the BGK scheme's performance to its 
derivation from a model to the collisional Boltzmann equa- 
tion. A tremendous advantage of a hydrocode designed on 
the basis of a solution to a model of such a fundamen- 
tal equation is that the scheme gets physical backing for 
its incorporation of viscous and heat conductive effects. 
Velocity moments of a time-dependent local solution to 
the BGK equation give hydrodynamical fluxes which carry 
both advective and dissipative terms - the BGK scheme 
does not decouple them into separate operations. Further- 
more the algorithm for solving the BGK equation in the 
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BGK scheme avoids the BGK equation's stiffness which 
would cause the time step of the scheme to be the colli- 
sion time of the gas instead of the much larger CFL time. 
Without such an algorithm to bypass the BGK equation's 
stiffness, the timestep would be prohibitively small for per- 
forming hydrodynamical simulations in real time. 

While it may seem complicated to solve for fluid evo- 
lution by following a changing distribution function (in- 
cluding the effects of collisions) in phase space, we believe 
that it is not. In spite of its looks, the BGK scheme is a 
computationally straightforward explicit scheme, and all 
parts of the algebra are in explicit closed form for a perfect 
gas. Even when a matrix is inverted for the connection be- 
tween the moments and parameters of a Maxwellian, the 
inverse is explicitly known. Unless readers were told they 
might not guess from inspection of the code that it origi- 
nated in gas-dynamic considerations. The code itself is just 
straightforward algebra. There are no iterative steps in the 
BGK scheme. There is no Riemann solver, either exact or 
approximate. No numerical sub-steps (e.g. Runge-Kutta 
steps) are used. 

Thus far we have not been concerned with the opti- 
mization of the code. Code development has focused on 
demonstrating that the code correctly simulates physical 
phenomena. Kim and Jameson (1998) compared the CPU 
time required just for the BGK scheme's flux computa- 
tion to the flux computation in two other high resolu- 
tion schemes which use flux-splitting methods: character- 
istic splitting using Roe averaging (Csplit) (Jameson 1996) 
and CUSP(Convective Upwind Split Pressure) splitting 
by Jameson (Kim and Jameson 1995). They find that the 
BGK scheme is less than twice as slow as the two other 
schemes. Note that in contrast to the flux computation in 
the other schemes, the BGK flux computation includes ex- 
tra arithmetical operations for Navier-Stokes terms. When 
the other schemes include Navier-Stokes terms, they are 
not spared this computational expense - it just appears 
outside of the flux computation. 
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Fig. 1. Early stage of time evolution of mass density profile for run with total energy conserving (ECS) scheme, 
including gravitational energy. Total time 201.19 (= 3.3 sound crossing times). 




Fig. 2. Early stage of time evolution of i-momentum density profile for run with ECS scheme. 
(=3.3 sound crossing times). 



Total time 201.19 




Fig. 4. Mass density profile for long run with ECS scheme after total time of 285739.54 (= 4706 sound crossing times). 




Fig. 5. Jf-momentum density profile for long run with ECS scheme after total time of 285739.54 (= 4706 sound 
crossing times). 
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Fig. 6. Lambda profile for long run with ECS scheme after total time of 285739.54 (= 4706 sound crossing times). 
Note that lambda becomes and remains constant after some time. This is entirely in accord with expectation of the 
long-term behaviour of a closed system: the temperature (i.e. j) reaches equilibrium and stays there. Compare this 
result with fig. |^ which shows long term heating of numerical origin using a different version of the code. 
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Fig. 7. Mass density profile for long run with energy source term (EST) scheme after total time of 261888.6 (= 4313 
sound crossing times). Note that the density profile has not settled to a constant profile with time. Compare with 
fig. [| which uses the ECS scheme. 



20 



A. Slyz & K.H. Prendergast: BGK with Gravity 




Fig. 8. X-momentum density profile for long run with EST scheme after total time of 261888.6 (= 4313 sound crossing 
times). 
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Fig. 9. Lambda profile for long run with EST scheme after total time of 261888.6 (= 4313 sound crossing times). 
Note that lambda descreases instead of remaining constant (as it should) at large time. The effect is due to numerical 
heating and is absent in the results from the totally conservative (ECS) scheme (cf. fig. @). 
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Fig. 10. Check to see how well final density profile from long run with ECS scheme matches the analytic solution for 
the equilibrium profile. The points should fall on a straight line if there is a perfect match to the analytic solution. 
See text (section |3.1|) for parameters of a least squares fit to the points. 
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Fig. 11. The initial conditions for the simulation of gas falling into a fixed external spherically symmetric gravitational 
potential well: uniform density (top left), lambda (top right), zero radial velocity (bottom left) and zero vertical velocity 
(bottom right). 




Fig. 12. The gravitational potential (top left), the radial force (top right) and the vertical force (bottom left), each 
as a function of w and z. They are fixed in time and space throughout the entire run. 
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Fig. 13. Results from simulations of the shock heating of gas falling into a fixed external spherically symmetric 
gravitational potential well. Time evolution of the density (top) and lambda (bottom) along two rays of the cylindrical 
grid. The cylindrical grid has 50 cells in w and 100 cells in z. 
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Fig. 15. Time evolution of the density (first column), lambda (second column) and velocity (third column) for the 
following sequence of times (in units of the initial free-fall time): 3.802D-03, 1.2671 , 2.723. 
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Fig. 16. Continuation of time evolution of the density (first column), lambda (second column) and velocity (third 
column) for the following sequence of times (in units of the initial free-fall time): 3.607, 4.8, 5.805. 



A. Slyz & K.H. Prendergast: BGK with Gravity 



29 




Fig. 17. Continuation of time evolution of the density (first column), lambda (second column) and velocity (third 
column) for the following sequence of times (in units of the initial free-fall time): 6.898, 8.067, 8.88. 
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Fig. 18. Lambda plotted once again at different z along w\ for i 9 — > *i 2 . This time the data for each of the four time 
instances is shown in a separate plot. The circles are the values of A in each cell and the thin lines are curves through 
these values. Note the sharpness of the shocks and the absence of near-shock oscillations. The shocks have not been 
sharpened in any way, nor have any oscillations been suppressed; all cells are shown in each plot. 
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Fig. 19. Density (top left), lambda (top right), radial velocity (bottom left) and vertical velocity (bottom right) after 
100,000 iterations which is equivalent to a total time in machine units of 36.3758. This amounts to about 560 free-fall 
times. Note that lambda is nearly constant and the velocity in both w and z is nearly zero, indicating that the code 
has (nearly) reached the expected equilibrium state. 



